source('solve.R')

funk1 <- function(method, dt)
{
	# R1: 0   -(b)-> T
	# R2: T   -(d)-> 0
	# R3: T+V -(e)-> Y+V
	# R4: Y   -(u)-> 0
	# R5: Y   -(p)-> Y+V
	# R6: V   -(c)-> 0

	# R4': Y  -(u)-> pV
	
	params <- list(
		#b=1, d=0.1, e=0.001, u=0.5, p=1000, c=0.5
		b=1, d=0.1, e=0.001, u=0.2, p=1000, c=10
	)
	X0 <- c(T=params$b/params$d, Y=10, V=10)

	# transition hazards
	transitions <- function(X)
		with(c(X, params), {
			c(b, d*T, e*T*V, u*Y, p*Y, c*V)
			#c(b, d*T, e*T*V, u*Y, c*V)
		})

	# stoichiometry matrix
	V <- with(params, {
		matrix(c(
			+1, 0, 0,
			-1, 0, 0,
			-1,+1, 0,
			# 0,-1, p,
			0,-1, 0,
			0, 0,+1,
			0, 0,-1
		), ncol=3, byrow=TRUE)
	})

	model.solve(100, X0, V, transitions, params, method, dt)
}

source('drawmodel.R')
drawmodel(
	function() funk1(model.euler, 0.01),
	20,
	function() funk1(model.tau_leaping, 0.01),
	c('funk1','model.euler','model.tau_leaping'),
	c('Target cells','Infected cells','Free virus')
)
